On an Optimal Number of Time Steps for a Sequential Solution of an Elliptic-Hyperbolic System

We consider a sequential approach for the solution of an elliptic-hyperbolic system of partial differential equations, which models a flow of two incompressible phases in porous media. The elliptic equation describes the pressure distribution in the domain, and the hyperbolic equation is the mass conservation equation for one of the phases. We propose to estimate an optimal number of the pressure updates using an analytical solution to a special 1D initial boundary value problem (IBVP) for the coupled system. We provide two procedures aimed at the estimation of an optimal set of time steps, and show that the resulting distribution of time steps yields better results than using equidistant time steps. We also show that the degree of coupling of the 1D IBVP can be quantitatively estimated using a normalized difference of the exact solution and its sequential approximation with a single time step. Introduction One well-established method for numerical solution of the equations for the multiphase flow in porous media is a sequential approach, which is based on a re-formulation of the mass conservation equations and Darcy’s law (see e.g. [1]). For the case of two incompressible phases with negligible capillary pressure and gravity effects (in this note we will restrict ourselves to this case), the re-formulated system of governing equations consists of the Laplace equation with variable coefficients for pressure, and the hyperbolic equation for the saturation of one of the phases. In the sequential approach, the solution to the coupled elliptic-hyperbolic system is sought in two stages: Firstly, the elliptic pressure equation with frozen coefficients is solved, and the corresponding Darcy velocity field is found; secondly, the hyperbolic saturation transport equation is solved over a certain period of time using a frozen Darcy velocity field, so that the coefficients of the elliptic equation can be computed; and the process repeats. A specific choice of the discretization technique for the elliptic and hyperbolic parts gives rise to a particular numerical method. Examples of the sequential method include IMPES-type methods [1, 2] and the method of streamlines [3]. One important question in the specification of the sequential approach is the choice of time instants at which the pressure field should be updated. In standard IMPES method, pressure is re-computed at each time step used for the solution to the hyperbolic transport equation, and the length of this common time step size is determined from CFL-like stability conditions [1, 2, 4]. In practice, this can require an excessive number of pressure updates for a simulation, which can make these methods computationally expensive. In the streamline method, one typically needs to estimate the number and length of pressure time steps before the actual computation. This decision is often based on engineering intuition. It has been noted (see e.g. [5]) that pressure is smoother in time than the saturation, and that one may use larger time steps for pressure updates than for the saturation transport. In [6], the authors obtained the convergence rates of a finite element method when the pressure time step is chosen as a fixed multiple of the saturation time step. Different pressure and saturation time steps are intrinsic to the method of streamlines, see e.g. [6] for an overview. An a posteriori CFL-like estimate for the length of pressure time steps in the streamline method is proposed in [7]. Finally, there is a variety of heuristic estimates for the length of pressure time steps available in the literature, based on maximum allowed pressure and saturation change, convergence criteria, etc.–see e.g. [1, 2]. To the best of our knowledge, currently there are no theoretical a priori estimates on the length of pressure time steps. This note is intended to fill this gap to some extent. To this goal, we restrict ourselves to a simple physical setting: we consider a flow of two incompressible fluids with negligible capillary pressure and gravity effects, whereas the relative permeability curves are straight lines. Moreover, let the flow be 1D, either in Cartesian geometry, or for cylindrically or spherically symmetric cases. For the resulting elliptic-hyperbolic system, we consider the following initial boundary-value problem (IBVP): piecewise constant initial data for saturation, and constant pressure boundary conditions. The initial and boundary conditions are chosen in such a way that a single shock-type solution is admissible. Then, we formulate an ordinary differential equation (ODE) for the motion of the interface, separating left and right states of constant saturation, and find its exact analytical solution. This is done in a way similar to Muskat [8]. For the case of Cartesian geometry, we present the dimensionless form of the IBVP of interest, and show that its single shock-type solution can be completely determined by specifying three dimensionless parameters: the mobility ratio, the initial interface position, and a parameter related to the shape of the Buckley–Leverett flux function. The sequential approach of the elliptic-hyperbolic system can be represented as a forward Euler integration of the ODE for the motion of the interface. Consequently, the sequential solution is a piecewise-linear approximation to the exact solution, and the vertices of this polygonal approximation correspond to the pressure time steps. We can estimate the accuracy of the polygonal approximation by e.g. considering its relative distance from the exact solution. Based on this error estimate, we propose two algorithms for an optimal time step selection: 1. Given an approximation error per time step which we agree to tolerate, we find the time instants at which the pressure update should be performed. 2. Given a total number of pressure updates, we find the time instants at which the pressure updates should be performed. These are chosen by minimizing the approximation error using a nonlinear optimization method–the downhill simplex method, see e.g. [9]. Finally, we propose to quantitatively estimate the degree of coupling of the 1D elliptichyperbolic IBVP using a normalized difference of the exact solution and its sequential approximation with a single time step. The functional dependence of the decoupling error with respect to the three dimensionless determining parameters is illustrated on a sample test problem. The rest of the paper is organized as follows. In Section 1, we formulate the mathematical problem of interest. In Section 2, we present the corresponding one-dimensional IBVP, integrate the ODE for the interface motion, and provide the admissibility conditions for the single shock-type solution, for Cartesian geometry and for cylindrically or spherically symmetric cases. For the case of Cartesian geometry, we present the dimensionless form of the IBVP and show that its solution is determined by three dimensionless parameters. In Section 3 we describe two algorithms for an optimal time step selection for the IBVP of interest. In Section 4 we define the decoupling error of the elliptic-hyperbolic system and illustrate this concept on a sample test problem. We end up with conclusions and outlook in Section 5. 1. Problem formulation In this work, we restrict ourselves to the case of two-phase incompressible fluids: oil and water. For the sake of simplicity we assume that both fluids differ in viscosity only and the rock is incompressible, and consider homogeneous media only, i.e. both porosity and permeability are constant. Furthermore, we neglect gravity and capillary forces. The mass conservation equations are